Empirical Functionals for Reduced Density Matrix Functional Theory 
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We present fully empirical exchange-correlation functionals to be used within reduced density 
matrix functional theory (RDMFT). These are of the popular J-K form, where the function of the 
occupation numbers that multiplies the Fock orbital term is written as a Fade approximant. The 
coefficients of the Fade are optimized for a test set of eight molecules, and then refined for a larger 
set of 35 molecules. Two different approaches were tried, either keeping the self-interaction terms, 
or by removing them explicitly from the functional. The functionals thus obtained involve very few 
parameters, but are able to outperform other RDMFT functionals, yielding correlation energies that 
are, on average, even slightly better than M0ller-Flesset MP2 theory. 
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Reduced matrix density functional theory (RDMFT) 
has been emerging as an excellent tool for the study of 
correlation in molecular systems. It is based on Gilbert's 
theorem, 1 that asserts a one-to-one correspondence be- 
tween the ground-state many-body wave function and 
the one-particle reduced density matrix, 7. Several the- 
oretical advantages are immediately evident from using 
7 instead of, e.g., the electronic density p as in standard 
density functional theory (DFT): i) both the kinetic en- 
ergy and the exchange energy can be trivially written as 
explicit functionals of 7; ii) consequently, the so-called 
correlation energy in RDMFT has a proper scaling with 
the strength of the Coulomb interaction, as it is free 
from contamination from the kinetic term. Therefore, 
one could expect that it is much easier to find good cor- 
relation functionals for RDMFT than to standard DFT. 

Indeed, the past years have witnessed the appearance 
of several such functionals. Most of them arc of the form 
(for closed-shell systems) 



3,k 



(1) 



i.e., they have the form of the usual Hartree-Fock (HF) 
exchange modified by the function f{nj, n^) of the occu- 
pation numbers. Functionals of the form of Eq. ([T]) are 
sometimes referred to as of J-K type, and involve only 
Coulomb (J) and exchange (K) type integrals over the 
natural orbitals. 

The first approximation to -Exc, introduced by Miiller 
in 1984,^ corresponds to the formula 

/M™-("j-,"fe) = NA^- (2) 
Later, Buijse and Baerends^ arrived at the same func- 



tional by modeling the exchange and correlation hole, 
while Goedecker and Umrigar (GU) considered a mod- 
ified version by explicitly removing the self-interaction 
(SI) terms. The extremely simple form of the Miiller 
functional yields the correct dissociation limit of dimers 
of open-shell atoms like II2 (where both DFT and HF 
fail), but overestimates substantially the correlation en- 
ergy^ of all the systems it has been applied to (including 
the IIEG)jii^i^ The GU form on the other hand, is more 
accurate in the calculation of the correlation energy j^i^ 
but fails dramatically at the dissociation limit .■AiS' 

Several other forms for f{nj,nk) exist in the market 
right now. The most precise for molecular systems seem 
to be the BBCS''^*' functional of Gritsenko, Pernal and 
Baerends, and the functionals of Piris.— These func- 
tionals give total correlation energies that are around 
17-20% from reference coupled-cluster CCSD(T) val- 
ues, just around a factor of two worse than M0ller- 
Plesset MP2 calculations Furthermore, atomization 
energies are basically of the same quality as DFT us- 
ing the B3LYP functional. Other properties of molecu- 
lar systems like ionization potentials j^^i^^'^^i^^ the chemi- 
cal hardnesSfi^iii dipole moments and static polarizabili- 
tiesr^i^i^ and vibrational frequencies were also calcu- 
lated with these functionals with very promising results. 

In this Article, we propose a fully empirical form 
— that we will denote by "ML" — for the function 
f(nj,nk). This approach, of using completely empiri- 
cal functionals, has an already long historjii^ in standard 
DFT, and yields some of the best functionals to datei^ 

We write f{nj,nk) as a general Pade approximant 



ao + aix + a2X + 



(3) 



where x — rijUk. We chose a Pade approximant as it is 
simple and one of the best choices for approximating a 



2 



rational function. Note that we multiplied the Pade by 
X to ensure that the contribution of completely empty 
states {uj = 0) to the exchange-correlation energy is zero. 
Furthermore, to recover the Hartree-Fock limit in the 
case of an idempotent density matrix, we force f^^{l) = 
1 by setting 



ao 



(4) 



4=1 



We also investigated the effects of removing the SI terms 
from the functionals (as suggested, e.g., in Refd), by 
trying the self-interaction corrected (SIC) functional 



/ 



ML-SIC 



S,,[rHn',)-n',] (5) 



No other constraints are imposed on the functional. This 
functional form was then optimized by a downhill simplex 
method, using as objective function 



S = 100 



1 ^ / 

i=l ^ 



/i;ML_ Rref 



(6) 



i.e., the average percentile error of the correlation energy 
for a test set of N molecules. For the reference correla- 
tion energy, i?™^ we used values computed with accurate 
coupled cluster theory [CCSD(T)] with the same basis 
set. 

The optimization of the Fade was performed in a very 
pragmatic way. First of all, we used a small number 
of parameters in the Fade. This reduces the problem 
of over-fitting of the functional, and simplifies the mini- 
mization procedure by keeping the dimensionality of the 
search space small. Furthermore, we used the Cartesian 
6-31G* Gaussian basis-set. This basis set is relatively 
small, allowing us to quickly perform the many calcula- 
tions needed to optimize the functional. The optimiza- 
tion was performed in two steps: i) The functional was 
minimized for a test set consisting on the eight of the 
smallest molecules of the G2 set)^ namely H2, hydro- 
gen fluoride, lithium hydride, water, ammonia, hydrogen 
chloride, hydrogen sulfide, and methane, ii) We then re- 
fined the Fade in the test set consisting of all closed-shell 
molecules in the G2-1 set (plus H2). This set includes 35 
molecules, with correlation energies ranging from around 
-0.02Hartree (LiH) to -O.SHartree (CO2). 

Our best functionals are summarized in Table ID Note 
that these functionals only depend on two parameters! 
In fact, we did not succeed in ameliorating the functional 
significantly by increasing their number. All the improve- 
ments were marginal (of less than 0.5%), and very often 
not consistent (improving the correlation energy of some 
subset of molecules but worsening others). We believe, in 
fact, that depending on just two parameters is a strength 
of our functionals, that in this way ally precision with 
simplicity. 



TABLE I: Parameters used to fit the functionals. Two ver- 
sions are given, one where the self-interaction terms are ex- 
plicitly excluded from the functionals, and another where they 
are not. Note that there are only two fitted parameters (ai 
and 61), as ao is obtained through Eq. Q. 



ao 



ai 



61 



ML 126.3101 2213.33 2338.64 

ML-SIC 1298.780 35114.4 36412.2 



TABLE 11: Error in the correlation energies for the closed- 
shell molecules of the G2 set calculated using the 6-31G* ba- 
sis set. The values in the second and third columns are in 
atomic units. See text for an explanation of the meaning of 
the columns. 



Method 



MuUer 
GU 
CGA 
BBCl 
BBC2 
BBC3 

PNOFO 
PNOF 
ML 

ML-SIC 



0.58 

0.28 

0.23 

0.31 

0.20 

0.074 

0.078 

0.12 

0.059 

0.062 



1.23 (C2CI4) 
0.79 (C2CI4) 
0.55 (C2CI4) 
0.75 (C2CI4) 
0.50 (C2CI4) 
0.27 (C2CI4) 
0.32 (SiCU) 
0.42 (SiCU) 
0.18 (pyridine) 
0.21 (pyridine) 



128.6% 
52.14% 
63.92% 
65.27% 
44.25% 
17.84% 
16.97% 
22.38% 
11.02% 
10.69% 



438.3% 
102.8% 
332.0% 
159.1% 
125.0% 
49.0% 
46.0% 
77.3% 
35.7% 
42.9% 



(Na2) 

(Na2) 

(Na2) 

(Na2) 

(Na2) 

(SiH2) 

(CI2) 

(F2) 

(Na2) 

(Na2) 



MP2 0.042 0.074 (C2CI4) 10.97% 35.7% (Lis) 



We then tested these functionals for all closed-shell 
molecules in the whole G2 set^ (119 molecules) using 
the same 6-31G* basis-set, and in the G2-1 basis set 
(35 molecules) using the better cc-pVDZ correlation con- 
sistent Dunning basis set. Results are summarized in 
Tables |TT] and IIIIl where we also included results ob- 
tained with M0ller-Plesser MP2 theory, and other of 
the most known RDMFT functionals, namely: Miiller— ; 
Goedecker and Umrigar (GU)^; Csanyi, Goedecker, and 
Arias (CGA)^^; Gritsenko, Pernal, and Baerends (BBCl, 
BBC2, and BBC3)^°; and the Piris functionals, both with 
the correction to avoid pinned states at one (PNOF) and 
without (PNOFO)." 

The meaning of the columns is the following 



A = 



An 
<5n 



1 ^ 

= max|£;c - E^"'^'> 

Er - E^°^ 

100 X max 



E^^ 



(7a) 

(7b) 
(7c) 



[The quantity S is defined as in Eq. ^.] 

The improvement of the functionals over the past years 
is quite remarkable. The best functionals that we have 
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TABLE III: Error in the correlation energies for the closed- 
shell molecules of the G2-1 set calculated using the cc-pVDZ 
basis set. The values in the second and third columns are in 
atomic units. See text for an explanation of the meaning of 
the columns. 



Method 


A 


^max 


5 




Miiller 


0.35 


0.57 (CI2) 


155% 


439% (Na2) 


GU 


0.13 


0.28 (CI2) 


46.4% 


120% (Na2) 


CGA 


0.16 


0.32 (S2) 


89.9% 


331% (Na2) 


BBCl 


0.19 


0.35 (CI2) 


72.8% 


181% (Na2) 


BBC2 


0.12 


0.24 (CI2) 


51.2% 


145% (Na2) 


BBC3 


0.048 


0.14 (CI2) 


21.4% 


68.9% (Na2) 


PNOFO 


0.045 


0.14 (CI2) 


18.0% 


44.6% (Na2) 


PNOF 


0.052 


0.16 (CI2) 


19.7% 


49.5% (CI2) 


ML 


0.026 


0.064 (N2) 


11.1% 


45.6% (Na2) 


ML-SIC 


0.023 


0.056 (N2) 


11.5% 


45.5% (Na2) 


MP2 


0.042 


0.074 (C2CI4) 


10.97% 


35.7% (Li2) 




FIG. 1: (Color online) Plot of the function f{nj,nk) for the 
ML, ML-SIC and the Miiller functional. The inset zooms on 
the particularly interesting region close to zero. 



available at the moment are literally an order of magni- 
tude better than the seminal Miiller functional, regard- 
less of the criterion used to judge them. We can also no- 
tice that our ML functionals are clearly the best RDMFT 
functionals. We especially remark the average percentile 
error S that shows an improvement of more than 60% 
with respect to the BBC3 result, being better than even 
MP 2 theory. The maximum errors are, unfortunately, 
still of a somewhat lower quality than MP2 theory. Note 
also that similar final results are obtained including or 
excluding the self-interaction terms in the functional. 

In Fig. [1] we plot the ML functionals compared to the 
simple Miiller form. It is clear that, in order to correct for 
the consistent over-correlation of the Miiller functional, 
our functionals are almost always smaller than the sim- 
ple square-root (with exception of a small region close 
to zero). We also note that the functionals never be- 
come negative (this was true for all forms that we tried, 
including those with more parameters). However, the 
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FIG. 2: (Color online) Histogram with the values of rijUk for 
all G2 closed-shell molecules. The values ofrijUk are obtained 
with the both with the ML (dashed red bars) and ML-SIC 
(empty green bars) functionals. Note that the left and right 
panels do not have the same horizontal or vertical scale. 



most striking fact is that the ML functionals are essen- 
tially linear except very close to zero, where their value 
drops quickly to zero (see inset of Fig. [T]). It is also im- 
portant to note that all functionals of reasonable quality 
that we found exhibited this kind of kink. 

To better understand our numerical results, we show, 
in Fig. [21 a histogram of the product rijUk for all G2 
closed-shell molecules, calculated both with the ML and 
ML-SIC functional. We can divide the product of the 
states in 4 different sets (we will use the nomenclature of 
Ref.lO): 

1. Products of weekly occupied orbitals (the large bar 
at zero). Note that in all cases, these weekly oc- 
cupied orbitals are not pinned at zero, but instead 
have a very small but finite occupation. 

2. Products of a weekly occupied orbital with a 
strongly occupied orbital. These are the prod- 
ucts whose values lie mainly between zero and 
~ 3 X 10"'^, but with a long tail that extends to 
UjUk = 0.25 (ML) and njUk = 0.15 (ML-SIC). 

3. Products of two strongly occupied orbitals. These 
have a tail that extends from njUk = 0.83 (ML) 
and HjUk = 0.95 (ML-SIC) to one. 

4. Products of two pinned states at one. For curiosity, 
we refer that both the ML and ML-SIC functionals 
yield exactly the same number of pinned orbitals 
for the molecules considered. 

The ML functional lead to occupation numbers (and cor- 
responding products rijUf^) with a much broader distri- 
bution than ML-SIC. This is not surprising, as the width 
of the distribution can be seen as a measure of "corre- 
lation". In the ML-SIC functional, a large contribution 
of the exchange-correlation energy (the self-interaction) 
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has been explicitly subtracted, so the resulting occupa- 
tion numbers will seem less correlated. We stress, how- 
ever, that from the point of view of the total correlation 
energy both ML and ML-SIC fare equally well, even if it 
apparently easier to construct SIC functionals. 

The analysis of the picture also makes clear a path to 
improve RDMFT functionals: to use different functional 
forms to approximate the exchange-correlation in each 
of the well-separated ranges corresponding to weakly- 
weakly, weakly-strongly, and strongly-strongly occupied 
states. The price to pay is a clear increase of the com- 
plexity of the functional. This is, in fact, the path already 
used by, for example, the BBC corrections.^ 

In conclusion, in this Article we proposed a very sim- 
ple empirical functional to be used within RDMFT. 
This functional is very precise, reaching in some respects 
the accuracy of more traditional quantum chemical ap- 
proaches, like MP2 theory. Perhaps surprisingly, we can 



reach the same level of average accuracy using functionals 
that are explicitly self-interacted corrected or not. This 
fact is further analyzed by looking at the basic ingre- 
dient of the functionals: the product njUk between the 
occupation numbers of two states. 
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